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Spin-dynamics techniques can now be used to study the deterministic time-dependent 
behavior of magnetic systems containing over 10 5 spins with quite good accuracy. This 
approach will be introduced, including the theoretical foundations of the methods of analy- 
sis. Then newly developed, improved techniques based upon Suzuki- Trotter decomposition 
methods will be described. The current "state-of-the-art" will be evaluated with specific 
examples drawn from data on simple magnetic models. The examination of dynamic criti- 
cal behavior will be highlighted but the extraction of information about excitations at low 
temperatures will be included. 

§1. Introduction 

The static behavior of physical systems near continuous phase transitions is 
characterized by a set of static critical exponents, which describe the critical be- 
havior of thermodynamic quantities such as the specific heat, the order parameter, 
the correlation length, and so on. One can thus define different universality classes, 
within which the critical exponents are identical. The numerical values of the critical 
exponents depend only on the symmetry of the order parameter, the dimensionality 
of the system and the range of interactions, but not on either the precise form of 
the model Hamiltonian or the lattice type. Likewise, the dynamic critical behav- 
ior is describable in terms of a dynamic critical exponent z, which depends on the 
conservation laws and which, in analogy to static critical phenomena, gives rise to 
different dynamic universality classes ^ . Our understanding of static critical behav- 
ior is now mature and has resulted largely from the investigation of simple model 
spin systems such as the Ising, the XY, and the Heisenberg model. These models 
are equally valuable for the investigation of dynamic critical behavior and dynamic 
scaling. Realistic models of magnetic materials can be constructed from these simple 
spin models; however, the theoretical analysis of experimentally accessible quantities, 
such as the dynamic structure factor, is usually too demanding for analytical meth- 
ods. Computer simulations are beginning to provide important information about 
dynamic critical behavior and material properties of model magnetic systems. 2 ) " 4 ) 
These simulations use model Hamiltonians with continuous degrees of freedom rep- 
resented by a three-component spin Sj with fixed length | Sj \ = 1 for each lattice site 
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j. A typical model Hamiltonian is then given by 

n = -jJ2 (s*sf + s]sf + \s*sf) -dJ2(s*)\ (i-i) 

<j,i> j 

where J is the coupling constant, < j,l > denotes a nearest-neighbor pair of spins, A 
is an anisotropy parameter, and D determines the strength of a single-site or crystal 
field anisotropy. (We use units in which Boltzmann's constant ks = 1.) 

The dynamics of the spins are governed by the coupled equations of motion 5 ) 

and the time dependence of each spin can be determined from the integration of these 
equations, where (hybrid) Monte-Carlo simulations of the model provides equilibrium 
configurations as initial conditions for Eq.(l-2). 

The most important quantity to be extracted from the numerical results is the 
dynamic structure factor 5(q, u) for momentum transfer q and frequency transfer 
to, which can be measured in neutron scattering experiments, and is given by 

/ + OO (If 
e ^t c k { ^ t) (1 . 3) 
n -°° V 2vr 

where R = rj — ri (rj and are lattice vectors), C fc (R, t) is the space- and 
time-displaced correlation function, with k = x,y, or z, defined as C k (R,t) = 

(Sj'm'm-is^msiHo)). 

Two practical limitations on spin-dynamics techniques are the finite lattice size 
and the finite evolution time. The finite time cutoff can introduce oscillations in 
S k (q,uj), which can be smoothed out by convoluting the spin-spin correlation func- 
tion with a resolution function in frequency, which is equivalent to the energy reso- 
lution in neutron-scattering experiments, yielding <S fc (q, uj). Finite-size scaling the- 
ory 3 ). 6 ) can De used to extract the dynamic critical exponent z: the divergence of 
the correlation length £ is limited by L and the dynamic finite-size relations are given 

by 

" S ~£?\ U) =GHuL',qL,6 u L>) (1-4) 

and 

^ = L-*W k (qL,8 u L*), (1-5) 

where x|(q) is the total integrated intensity and cZ^ is a characteristic frequency, 
defined as 

f^m - k diU 1 u 

/ S L \c { ,u J )— = -x k M). (1-6) 

To speed up the numerical integration of Eq.(l-2) it is desirable to use the 
largest possible time step; however, with standard methods the size of the time 
step is severely limited by the accuracy within which the conservation laws of the 
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dynamics are obeyed. It is evident from Eq.(l-2) that the total energy is conserved, 
and if, for example, D = and A = 1 (isotropic Heisenberg model) the magnetization 
M = J2j Sj is also conserved. For the anisotropic Heisenberg model, i.e., A / 1 or 
D ^ only M z is conserved. Conservation of spin length and energy is particularly 
crucial, because the condition |Sj| = 1 is a major part of the definition of the model 
and the energy of a configuration determines its statistical weight. It would therefore 
also be desirable to devise an algorithm which conserves these two quantities exactly. 

The remaining sections of this paper are organized as follows. In Section 2 
we describe integration methods, including a newly developed technique based on 
Suzuki- Trotter decompositions of exponential operators. 7 )' 8 ) In Section 3 we discuss 
two examples of physical systems, namely the two-dimensional XY model and the 
three-dimensional Heisenberg antiferromagnet, and comparison with experiment and 
theory. The purpose of these examples is to show just how far the state-of-the-art 
has developed in producing and analyzing spin-dynamics data. In Section 4 we give 
a brief summary of this paper. 

§2. Integration methods 

2.1. Predictor- corrector method 

Predictor-corrector methods have been quite effective for the numerical integra- 
tion of spin equations of motion; however, in order to limit truncation errors small 
time steps St must be used with at least a 4th-order scheme. The predictor step of 
the scheme used here is the explicit Adams-Bashforth four-step method 9 ) and the 
corrector step consists of typically one iteration of the implicit Adams- Moulton three- 
step method. 9 ) This predictor-corrector method is very general and is independent 
of the special structure of the equations of motion (see Eq.(l-2)). The conservation 
laws discussed earlier will only be observed within the accuracy set by the truncation 
error of the method. In practice, this limits the time step to typically A > St = 0.01/ J 
in d = 3, where the total integration time is typically 600/ J or less. 

2.2. Suzuki- Trotter decomposition methods 

The motion of a spin may be viewed as a precession of the spin S around an 
effective axis ft which is itself time dependent. The lattice can be decomposed into 
two sublattices such that a spin on one sublattice precesses in a local field ft of 
neighbor spins which are all located on the other sublattice. For the Hamiltonian in 
Eq.(l-l) there are only two such sublattices if the underlying lattice is simple cubic. 

To illustrate the method, we consider first the D = case. The basic idea of the 
algorithm is to rotate a spin about its local field ft by an angle a = \ ft\5t, rather than 
directly integrate Eq.(l-2). This procedure guarantees the conservation of the spin 
length | S | and energy to within machine accuracy. Denoting the two sublattices by A 
and B, respectively, we can express the local fields acting on the spins on sublattice A 
and B as 1?^[{S}] and l?g[{S}], respectively. In a more symbolic way, we denote y as 
a complete spin configuration, which is decomposed into two sublattice components 
y_4 and yg, i.e. y = {yAiVs)-, and we denote by matrices A and B the generators 
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of the rotation of the spin configuration yj± on sublattice A at fixed yg and of the 
spin configuration on sublattice B at fixed y^, respectively. The update of the 
configuration y from time t to t+St is then given by an exponential (matrix) operator 

y{t + 5t) = e^ A+B ^y(t). (2-1) 

The exponential operator in Eq.(2T) rotates each spin of the configuration and it 
has no simple explicit form, because the rotation axis for each spin depends on 
the configuration itself; however, the set of equations of motion for spins on one 
sublattice reduces to a linear system of differential equations if the spins on the 
other sublattice are kept fixed and the operators e ASt and e BSt which rotate yjs, at 
fixed ys and ys at fixed y^, respectively, do have a simple explicit form. 8 ) Thus an 
alternating update scheme may be used, i.e., we rotate y^ at fixed yjg and vice- versa. 
The alternating update scheme amounts to the replacement e( A+B ) &t — > e A5t e BSt j n 
Eq.(2-1), which is only correct 10 ) up to order (5t) 2 . The magnetization will therefore 
only be conserved up to terms of the order 5t (global truncation error). To decrease 
truncation error and thus to improve the conservation, one can employ mth-order 
Suzuki- Trotter decompositions of the exponential operator in Eq.(2-1), namely 10 ) 

u 

e (A+B)8t = -Q e piA8t/2 ePi B8t ePi A5t/2 + Q^t m+1 ) (2-2) 
i=l 

where u = 1 for 2nd-order, u = 5 for 4th-order, and u = 15 for 8th-order and the 
parameters Pi are given by Suzuki and Umeno. 10 ) 

The additional computational effort needed to evaluate higher order expressions 
can be compensated to some extent by using larger time steps. The inclusion of 
next-nearest neighbor bilinear interactions on a simple cubic lattice can be treated 
within the above framework if the lattice is decomposed into four sublattices. 

This approach can also be extended to the case but in contrast to the 

isotropic case, the equation of motion for each individual spin on each sublattice is 
nonlinear. In practice, the best form of solution is via iterative numerical methods. 
In order to perform a rotation operation in analogy to the isotropic case we identify an 
effective rotation axis Qj = Qj — D ^0, 0, Sj(t) + Sj(t + 5t)^j , such that the condition 

for energy conservation is rewritten in the form fij ■ (Sj(t + 5t) — Sj(t)) = 0. Since 
the rotation requires knowledge of Sj at the future time t + 8t, we use an iterative 
procedure starting from the initial value Sj(t + 5t) = Sj(t) + (Qj x Sj(t)) z 5t and 
performing several updates according to the decompositions given by Eq.(2-2) in 
order to improve energy conservation. Both the degree of conservation and the 
execution time depend to some extent on the number of iterations used. 

For a quantitative analysis of the integration methods outlined above we restrict 
ourselves to the Hamiltonian given by Eq.(lT) for A = 1 in d = 3 and J > 0. 
The underlying lattice is simple cubic with L = 10 lattice sites in each direction and 
periodic boundary conditions. In order to compare the different integration methods 
we investigate the accuracy within which the conservation laws are fulfilled. The 
initial configuration is a well equilibrated one from a Monte-Carlo simulation for 
A = 1 at a temperature T = 0.8T C for D = and D = J, where T c refers to the 
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critical temperature of the isotropic model (D = 0). Fig.l shows the magnetization 
conservation for the 4th- and 8th-order decomposition methods, both with St = 
0.1/ J, for D = and Fig. 2 shows the energy conservation for different methods for 
D = J. In the latter case the iterative nature of all four methods gives rise to a 
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Fig. 1. Magnetization m(t) = |M(£)|/L 3 per spin for different order decomposition schemes for 
D — and time step St — 0.1/. J: (solid line) 4th-order scheme; (dashed line) 8th-order method. 
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Fig. 2. Energy e(t) = E(t)/(JL 3 ) per spin for different order decomposition schemes for D — J: 
(solid line) predictor-corrector method; (dot-dashed line) 2nd-order scheme; (dashed line) 4th- 
order scheme; (dotted line) 8th-order method. The number of iterations performed are marked 
next to each line. 



basically linear energy change. A single integration step using the 2nd-, the 4th- and 
the 8th-order scheme is respectively about 2 times faster, 2.5 and 9 times slower than 
the predictor-corrector method; the speedup of the decomposition methods comes 
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from the much larger St that can be used. For D = it is still feasible to use the 
4th-order method with St = 0.2/ J, which corresponds to an eightfold speedup as 
compared to the predictor-corrector method. The 8th-order method improves the 
conservation significantly but at the cost of greatly increased execution time. 

§3. Examples of physical systems 

3.1. Two-dimensional XY model 

The two-dimensional XY model can be described by the Hamiltonian in Eq.(lT) 
with A = and D = 0. At the critical temperature 2 ) Tkt = 0.700(5) J, the model 
undergoes an unusual phase transition to a state with bound, topological excitations 
(vortex pairs), and the static properties are consistent with the predictions of the 
Kosterlitz-Thouless theory. n ) For T < Tkt the model is critical, i.e. the correlation 
length is infinite, but there is no long-range order, and the spin-spin correlation 
function decays algebraically with distance, with an exponent rj that varies with 
temperature. 

In our studies of the dynamics of this model 2 ) , we used L x L lattices with peri- 
odic boundary conditions for 16 < L < 192 and several values of T. The equations of 
motion (see Eq.(l-2)) were integrated using a 4th-order predictor-corrector method, 
with a time step of St = 0.01/ J, to a maximum time t max = 400/ J. Between 500 
and 1200 equilibrium configurations were used for each lattice size and temperature. 
We were limited to the [10] reciprocal lattice direction, i.e. q = (q, 0) and (0, q). 

For T < Tkt, the in-plane component S x (q,ui) exhibits very strong and sharp 
spin-wave peaks. As T increases, they widen slightly and move to lower lu, but re- 
main pronounced even just above Tkt- For increasing momentum they broaden and 
rapidly lose intensity. Well above Tkt, the spin- wave peak disappears in S x (q, u), as 
expected, and we observe a large central peak instead. Besides the spin-wave peak, 
S x (q,uj) exhibits a rich low-intensity structure, which we interpret as two-spin- wave 
processes (see Fig. 3). Furthermore, S x (q,uj) shows a clear central peak, even below 
Tkt, which becomes very pronounced towards Tkt- Neither this strong central peak 
nor the additional structure are predicted by existing analytical calculations. The 
out-of-plane component S z (q,uj) is much weaker than S x (q,co), except for large q. 
The very sharp spin-wave peaks at low temperatures allowed us to determine the 
dispersion curves with great accuracy. Our estimated value of the dynamic criti- 
cal exponent is z = 1.00(4), in agreement with the theoretical prediction of z = 1. 
The line shape of S x (q,Lo) is not well described by either Villain's 12 ) or Nelson and 
Fisher's 13 ) prediction (the latter agrees qualitatively with our data only for large q) 
(see Fig.4). Moreover, these predictions do not describe the additional structure in 
S x (q,uj), including the central peak. 

3.2. Three-dimensional Heisenberg antiferromagnet and RbMnF^ 

RbMnF3 is a good physical realization of an isotropic three-dimensional Heisen- 
berg antiferromagnet, described by the Hamiltonian in Eq.(lT) with A = 1, D = 
and J < 0. Early experimental studies [see references in Ref.14)] showed that the 
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Fig. 3. Low-intensity structure in S x (q,u>) for T=0.6J, L=192 and q = 7r/32. Vertical arrows show 
the location of two-spin- wave peaks formed by spin waves of small momentum q < 4(2tt/L). 
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Fig. 4. Comparison of the line shape of S x (q,uS) with theoretical predictions. Data are at T = Tkt, 
L = 128, and q — tv/32 (thick line), normalized as cqS x (q,u) /\ x (q), where c is the spin-wave 
velocity. The two thin lines represent the predictions by Nelson and Fisher 13 ' (continuous line) 
and by Villain 12 ' (dashed line), both with rj = 0.25. The inset shows the data and predictions 
on a log- log plot that includes large values of ui. 



Mn 2+ ions, with spin S = 5/2, form a simple cubic lattice structure with a nearest- 
neighbor exchange constant \ J exp \ = 0.58(6) meV and a second-neighbor interaction 
constant of less than 0.04 meV, both defined using the normalization as in Eq.(l-l). 
Magnetic ordering with antiferromagnetic alignment of spins occurs below the criti- 
cal temperature T c = 83K. The magnetic anisotropy is very low, about 6 x 10~ 6 of 
the exchange field, and no deviation from cubic symmetry was seen at T c . 

In our simulations 14 ^, we used simple cubic lattices with 12<L<60atT< 
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T c = 1.442929(77)| J\. Numerical integrations of the coupled equations of motion 
were performed to a maximum time t max < 1000] J] -1 , using the algorithm based on 
4th-order Suzuki- Trotter decompositions of exponential operators, with a time step 
St = 0.2 1 J| _1 • As many as 7000 initial configurations were used, although for large 
lattices this was reduced to as few as 400. We were limited to the [100], [110] and 
[111] directions, i.e. q = (q, 0, 0), (q,q,0), (q,q,q) and the equivalent momenta. For 
this model the order parameter is not conserved and the dynamic structure factor 
cannot be separated into a longitudinal and a transverse component. Henceforth we 
will use the term dynamic structure factor S(q,to) and characteristic frequency Q m 
to refer to the average. 

We compare our results with the recent neutron scattering data of Coldea et 
al. 15 ) Although the compound RbMnF3 is a quantum system, and our simulations 
are for a classical Hamiltonian, it has been shown that quantum Heisenberg sys- 
tems with large spin values (S > 2) behave as classical Heisenberg systems where 
the spins are vectors of magnitude \J S{S + 1) with the same interaction strength. 
Since our classical spins were vectors of unit length, to preserve the Hamiltonian 
and the dynamics, a normalization of the interaction strength and the frequency 
transfer are needed, i.e. J = J exp S{S + 1) and uj exp = J exp ^S{S + 1) w/J. For our 
comparison, the experimental T- and w-dependent population factor was removed 
from the experimental data, and the lineshapes from our simulation were convoluted 
with the experimental Gaussian resolution function in energy. Fig. 5 shows the di- 
rect comparison of S(q,u>) in the [111] direction from simulation and experiment 
for q = 2-7r(0.08) [in this notation the Brillouin zone edge in the [111] direction cor- 
responds to q = 27r(0.25)], at T = 0.894T C and at T c . Below T c , renormalization 
group theory 16 ^ (RNG) predicts a spin-wave peak and a central peak in the longi- 
tudinal component of S(q,uj); however, at T c , both RNG 17 ) and mode coupling 18 ) 
(MC) theory predict only the presence of a spin-wave peak, while the experiment 
and the simulation find a spin-wave peak and a central peak at T = T c as well (see 
Fig. 5(b)). At low temperatures the central peak has very low intensity and the dom- 
inant structures are very narrow and sharp spin-wave peaks, from which accurate 
dispersion curves could be found. The dispersion curve for small q changes from a 
linear behavior at low T to a power-law relation as T — ► T c . 

The dynamic critical exponent z was extracted from the slope of a log(<D m ) vs 
log(L) plot (see Eq.(l-5)) corresponding to the [100] direction, using lattices in the 
asymptotic-size regime (L > 30), and keeping qL and b^L z constant. For 5^ = 0, 
we find z = 1.45(1) for n = qL/(2n) = 1 and z = 1.42(1) for n = 2. Using 
b~u) 7^ requires an iterative procedure and the converged values that we obtained 
are z = 1.43(1) for n = 1 and z = 1.42(1) for n = 2. Hence, our final estimate is 
z = 1.43(3), which is in agreement with the recent experimental value z = 1.43(4), 
and slightly lower than the theoretical prediction z = d/2 = 1.5. 

§4. Summary 

We have shown how spin-dynamics techniques can be used to study critical and 
low-temperature magnetic excitations using simple classical spin models that have 
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Fig. 5. Comparisons of lineshapes obtained from fits to simulation data for L = 60 (solid line) and 
the experiment (open circles) for q — 27r(0.08) in the [111] direction, at (a) T = 0.894T C and (b) 
T = T c . The dot-dashed line in (b) is a fit to the experimental data which is compared to the 
predictions of renormalization group (RNG) and mode coupling (MC) theory. The horizontal 
line segment in each graph represents the 0.25meV resolution in energy (full-width at half- 
maximum) . 



true dynamics, governed by equations of motion. The solution of these equations 
is generally possible through the use of algorithms based on Suzuki- Trotter decom- 
positions of exponential operators and we compare their relative performance with 
each other and with a predictor-corrector method. As examples of interesting phys- 
ical systems, we studied the two-dimensional XY model and the three-dimensional 
Heisenberg model. We determined dynamic structure factors and through a finite- 
size scaling we estimated the dynamic critical exponent of these systems. We have 
also made comparisons with theoretical predictions and experimental data. 
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Fig. 6. Finite-size scaling plot for u> m (with qL =const, <LL Z =const) for the analysis with and 
without a resolution function. For the former case, the data used correspond to the converged 
values of z, for n = 1, 2. The error bars were smaller than the symbol sizes. 
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